In recent years, there has been a resurgence in developing new solver technologies for addressing highly complex and large-scale flow simulations on specialized parallel and multicore architectures in a very effective manner. Methods such as algebraic multigrid (AMG), Krylov recycling (e.g., deflation, Krylov-secant) and extensions to two-stage preconditioners (e.g., GPR) have been incorporated into the latest solver advances.
Jenny, et al., U.S. Pat. No. 6,823,297, teaches a multi-scale finite-volume (MSFV) approach for solving elliptic or parabolic problems such as those found in subsurface flow simulators. Wallis and Tchelepi, US2007010979, use an algebraic cascading class to reduce cell nodes to two or more domains and graphing the model. Lee, et al., US2008208539, utilize a multi-scale finite-volume (MSFV) method to deal with nonlinear immiscible three phase compressible flow in the presence of gravity and capillary forces. Poulisse, et al., WO2008006851, use algebraic methods to generate a polynomial ring that describes the polynomial relationship between the variables. Finally, Carruthers and Neufeld, US2008262809, model the migration of a reactant in a petroleum reservoir using a mesh comprising a plurality of nodes to calculate one or more variables in the model.
A recent survey of methods to perform a static characterization of reservoirs is given by Hovadik and Larue (2007). Percolation algorithms offer efficient ways to perform qualitative predictions on random media (Torquato, 2001; Bollobas and Riordan, 2006). In particular, the Hoshen-Kopelman algorithm (1976) is one of the most popular approaches to find critical paths, percolation thresholds, and conductive clusters. In many pore network studies, the random media is represented by a set of sites (nodes) and bonds (links) labeled according to the presence or absence of a particular property (e.g., conductance, reaction condition). The determination of the number, size and shape of clusters in a percolation graph is critical to perform reliable statistics on the conductive pattern in such random systems. Most of the clustering analysis is numerically estimated with the aid of the H-K algorithm on lattice environments (Hoshen and Kopelman, 1976). However, an efficient version of this clustering algorithm was recently proposed in (Al Futaisi and Patzek, 2003) that is suitable for arbitrary geometries.
Nodes can be visited in an either depth-first search or breadth-first search fashion (Jungnickel, 2008) giving rise to an algorithm of linear complexity in terms of the sum of the number of nodes and edges. Moreover, this algorithm can be implemented with slight modifications to the better known family of spanning trees algorithms (Babalievski, 1997). Aggregation or agglomeration can be used to gather strongly connected coefficients into a single set (aggregate). The reduction of these aggregates leads to the construction of coarse operators that are suitable to reduce the computational burden of the original problem (see e.g. Fish and Belsky 1997; Brezina et al., 2005; Muresan and Notay 2008). Aggregation is usually achieved by mapping the coefficients into a graph and then finding a maximal independent set consisting of no adjacent nodes—or equivalently, a collection of edges that do not share a common node.
High coefficient contrasts commonly arise, for example, in multilayered geological formations composed of different type of rocks (e.g., sandstones, shales, carbonates, limestones) (Torquato, 2001). Understanding the role of these systems characterized by highly heterogeneous media has been the subject of intensive research to account for all possible scales present in the problem (see e.g., Hou et al., 1999; Lackner and Menikoff, 2000; Tchelepi et al., 2005; MacLachlan and Moulton, 2006). Nevertheless, many of these approaches disregard the fact that the random media often contain connected regions defining preferential solution paths (e.g., fluid flow, pressures, electrical conductance, resistivity or transmissibility depending on the application) that may be helpful to improve the solution process. From the linear algebra standpoint, some authors have been particularly interested in establishing a “downwind” or “crosswind” criteria to adapt solver algorithms to follow directions with strong coupling (Wang and Xu, 1999; Kim et al., 2003). In the general case of pressure systems arising from IMPES formulations, a symmetric formulation may be preferred (see e.g., Naben and Vuik, 2006).M=(I−AM1)M2(I−M1A)+M1 
A preprocessing scaling step before applying the whole main preconditioner step has been proposed but has not been shown in use with a preconditioner system (Graham and Hagger, 1999; Aksoylu and Klie, 2008). A matrix compensation algorithm (Lu, 1995; Huang et al. 2007) can be optionally applied to discard weakly connected coefficients. The homogeneous system can be described by relaxation sweeps; (see e.g., Brezina et al., 2005; Huang et al. 2007). Spectral preconditioners (Giraud and Gratton, 2006), can be seen as a powerful complement to improve the performance of existing preconditioners; (see e.g. Frank and Vuik, 2001; Naben and Vuik, 2006; Aksoylu and Klie, 2008)
Multilevel algebraic procedures cannot address assumptions present in developing reservoirs. These traditional algebraic methods do not accurately reflect positive definiteness, diagonal dominance, and the symmetry of linear systems operators found in subterranean formations and production data. Although current algebraic systems can be used to efficiently describe moderately sized and relatively well posed problems, the cannot be applied to highly heterogeneous and anisotropic media for large scale petroleum reservoir simulations.